!---------------------------GRELL CUMULUS SCHEME---------------------------
subroutine CUPARTH(mynum, mgmxp, mgmyp, mgmzp, m1, m2, m3,        &
     ia, iz, ja, jz, i0, j0,                                      &
     maxiens, iens,ngrid, ngrids_cp,                              &
     DTIME, time, UA, VA, WA, THETA, PP, PI0, DN0, RV, HT, RTGT,  &
     THT, RTT, PT, OUTTEM, OUTRT, PRECIP, sgrell,                  &
     !Variables to save cumulus properties
     ierr4d, jmin4d, kdet4d, k224d, kbcon4d, ktop4d, kpbl4d,      &
     kstabi4d, kstabm4d, xmb4d, edt4d, zcup5d, pcup5d, enup5d,    &
     endn5d, deup5d, dedn5d, zup5d, zdn5d,  prup5d, &	!57 Lufla ???
     clwup5d, &	!58 Lufla ???
     tup5d,                 &
     !Variaveis to save in CUP_DIRECTION2
     upmf, dnmf, xierr, xktop, xkbcon,xk22, xjmin, xkdt, xiact_p, xiact_c)

  !srf   mgmxp, mgmyp, mgmzp sao usadas alocar memoria para as
  !      variaveis da parametrizacao do Grell.

  ! USE Modules for Grell Parameterization
  use mem_grell_param, only : maxens,  & !INTENT(IN)
       maxens2,                        & !INTENT(IN)
       maxens3,                        & !INTENT(IN)
       ensdim,                         & !INTENT(IN)
       icoic                             !INTENT(IN)

  use mem_scratch2_grell

  implicit none
  integer :: mgmxp,mgmyp,mgmzp,ngrid,ngrids_cp
  integer :: maxiens !,maxens,maxens2,maxens3,ensdim
  integer :: iens   ! Local Variable

  !  this is the lowest level, then comes ensemble 1!
  !       ensdim=maxiens*maxens*maxens2*maxens3 !Ensemble dimension
  !INTEGER icoic
  !
  !srf- ICOIC is used for choice a specific closure
  ! icoic = 0 -> ensemble (all closures)  [en]
  ! icoic = 1 -> Grell                    [gr]
  ! icoic = 4 -> low level omega          [lo]
  ! icoic = 7 -> moisture convergence     [mc]
  ! icoic =10 -> like Fritsch Chappel or Kain Fritsch [sc]
  ! icoic =13 -> Arakawa-Schubert         [as]
  !PARAMETER (icoic=1)

  ! Defined in mem_grell_param

!!!!!!!!!!!!

  !-------------
  ! Constantes:
  real, parameter :: rgas = 287., cp = 1004., rm = 461., p00 = 1.e5,  &
       tcrit = 273.15, g = 9.80, cpor = cp/rgas, PKDCUT = 75.

  !------------------------------------ RAMS vectors:
  integer m1,m2,m3,ia,iz,ja,jz,i0,j0,ibcon,j1,j2,mynum
  real time
  real, dimension(m1,m2,m3) :: ua, va, wa, theta, pp, pi0, dn0,   &
       rv, tht, rtt, pt, OUTTEM, OUTRT
  REAL, DIMENSION(m1,m2,m3) :: SGRELL
  real, dimension(m2,m3)    :: ht,rtgt,PRECIP

  real, dimension(m2,m3)    :: upmf, dnmf, xierr, xktop, xkbcon, xjmin,  &
       xkdt, xiact_p, xiact_c, xk22

  !-------Salva parametros da CUP para uso no transporte convectivo:
  integer, dimension(mgmxp,mgmyp,maxiens,ngrids_cp)    :: ierr4d, jmin4d,  &
       kdet4d, k224d, kbcon4d, ktop4d, kpbl4d, kstabi4d, kstabm4d

  real, dimension(mgmxp,mgmyp,maxiens,ngrids_cp)       :: xmb4d, edt4d

  real,dimension(mgmzp,mgmxp,mgmyp,maxiens,ngrids_cp) :: enup5d,  &
       endn5d, deup5d, dedn5d, zup5d, zdn5d, &   !p_lw5d, &
       prup5d, &  !Lufla
       clwup5d, & !Lufla
       tup5d, &   !Lufla
       zcup5d, &  !Lufla
       pcup5d

  !------------------------------------- variaveis locais:

  integer :: kk, istart, iend, i, j, k, mix, mjx, mkx, kr, m, ipr, jpr
  real    :: vspeed, dp, dtime, dq, cpdTdt, exner

  ! Init. arrays P and PO with zeros - ALF

  do k=1,mgmzp
     do i=1,mgmxp
        P(i,k)  = 0.
        PO(i,k) = 0.
     enddo
  enddo

  !----------------------------------------------------------------------
  ISTART = ia
  IEND   = iz
  j1     = ja
  j2     = jz
  MKX    = m1 - 1    !MKX nao deve ser igual a m1
  MIX    = m2            
  MJX    = m3            

  !---- Coordenadas para escrita ascii
  !      ipr=19 - i0
  !      jpr=19 - j0
  ipr=0 - i0
  jpr=0 - j0

  do j=1,m3  ! loop em todo dominio para passar informacoes da fronteira
     do i=1,m2 ! dos nodes
        massflx(i,j)     = dnmf(i,j)
        iact_gr(i,j)     = int(xiact_c(i,j))
        iact_old_gr(i,j) = int(xiact_p(i,j))
        !        print*,massflx(i,j),iact_gr(i,j),iact_old_gr(i,j)
     enddo
  enddo

  ! Loop externo : j

  do J=j1,j2

     !DO 20 I = ISTART,IEND
     do I = ISTART,IEND
        AA0(I)=0.
        XLAND(I,J) = 0. ! land/water flag - not in use
!20      CONTINUE
     enddo

     do I = ISTART,IEND
        IACT_GR(I,J)     = 0    !verificar se isto esta' correto
        IACT_OLD_GR(I,J) = 0    !verificar se isto esta' correto
        KDT(I,J)         = 0
        PRECIP(I,J)      = 0.
        MASSFLX(I,J)     = 0.   !verificar se isto esta' correto
     enddo

     !--- Prepare input, erase output

     do I = ISTART,IEND        
        KDET(I)  =2
        PRET(I)  =0.
        MCONV(I) =0.
        UMEAN(I) =0.
        VMEAN(I) =0.
        PMEAN(I) =0.
     enddo

     !------- Transfere valores do RAMS para o eschema
     do K=1,MKX
        kr = K + 1          ! nivel K da grade do Grell corresponde ao
        ! nivel K + 1 do RAMS
        do I = ISTART,IEND

           TER11(I)= ht(i,j)
           PSUR(I) = .5*( ((pp(1,i,j)+pi0(1,i,j))/cp)**cpor*p00 +  &
                ((pp(2,i,j)+pi0(2,i,j))/cp)**cpor*p00 )*1.e-2
           ! Pressure in mbar

           PO(I,K) = ((pp(kr,i,j)+pi0(kr,i,j))/cp)**cpor*p00*1.e-2
           ! Pressure in mbar
           US_Grell(I,K) = .5*( ua(kr,i,j) + ua(kr,i-1,j) )
           VS_Grell(I,K) = .5*( va(kr,i,j) + va(kr,i,j-1) )
           OMEG(I,K)   = -g*dn0(kr,i,j)*.5*( wa(kr,i,j)+wa(kr-1,i,j) )

           T(I,K)  = theta(kr,i,j)*(pp(kr,i,j)+pi0(kr,i,j))/cp
           Q(I,K)  = rv(kr,i,j)

           !        Calcula tendencia projetada na temperatura em funcao 
           !        das tendencias de theta e PI : cp*T=Pi*Theta
           exner= pp(kr,i,j)+pi0(kr,i,j)

           !        cpdTdt= exner*tht(kr,i,j) + theta(kr,i,j)*pt(kr,i,j)
           cpdTdt  = exner*tht(kr,i,j)
           ! assumindo PT(KR,I,J) << exner*THT(KR,I,J)/theta

           !        Temperatura projetada se a conveccao nao ocorrer
           TN(I,K) = T(I,K) + ( cpdTdt/cp )*dtime

           !        Umidade projetada se a conveccao nao ocorrer
           QO(I,K) = Q(I,K) +   rtt(kr,i,j)*dtime

           !srf-----print-------
           if(j.eq.jpr   .and. i.eq.ipr) then
              if(k.eq.1) then
                 write(2,'(a1,78a1)') ' ',('-',m=1,78)

                 write(2,1111) IPR,JPR,MKX,ter11(i),time/3600.,psur(i),dtime
1111             format(1x,3i4,4E12.5)
              endif
              write(2,1112) PO(i,k), T(I,K), Q(I,K), TN(I,K), QO(I,K),  &
                   US_Grell(i,k), VS_Grell(I,k), omeg(i,k)
1112          format(1x,9E12.5)
           endif

           !         if(j.eq.jpr   .and. i.eq.ipr) then
           !	  if(k.eq.1) then
           !	   print*,'I K PO T Q theta rv dn0'
           !	  endif
           !          write(6,'(2i4,6F12.4)') i,k,PO(i,k),T(I,K),Q(I,K),
           !    &                   theta(kr,i,j),rv(kr,i,j),dn0(kr,i,j)
           !         write(6,'(2i4,5F12.4)') i,k,TN(I,K),QO(I,K),cpdTdt/cp
           !     &    ,rtt(kr,i,j),dtime
           !          if(k.eq.mkx) write(6,'(a1,78a1)') ' ',('-',m=1,78)
           !         endif
           !srf-----print-------

           !------- Atribuicoes do eschema

           P(I,K)  = PO(I,K)
           !srf	 PSUR(I) = 0.5*(PO(I,1)+PO(I,2))
           if((PSUR(I)-P(I,K)).gt.150.and.P(I,K).gt.300.)then
              DP       = -.5*(P(I,K+1)-P(I,K-1))
              UMEAN(I) = UMEAN(I)+US_Grell(I,K)*DP
              VMEAN(I) = VMEAN(I)+VS_Grell(I,K)*DP
              PMEAN(I) = PMEAN(I)+DP
           endif

           if(TN(I,K).lt.200.)    TN(I,K) = T(I,K)
           if(QO(I,K).lt.1.E-08)  QO(I,K) = 1.E-08

           OUTT(I,K) = 0.
           !Tendencia no campo de temperatura associada aos cumulus
           OUTQ(I,K) = 0.
           !Tendencia na razao de mist. de vapor d'agua assoc. aos cumulus
           OUTQC(I,K) = 0.
           !Tendencia na razao de mistura de agua de nuvem e/ou gelo
           ! associada aos cumulus
        enddo
     enddo

     do I = ISTART,IEND
        UMEAN(I)=UMEAN(I)/PMEAN(I)
        VMEAN(I)=VMEAN(I)/PMEAN(I)
        VSPEED=sqrt(UMEAN(I)*UMEAN(I)+VMEAN(I)*VMEAN(I))
        DIRECTION(I)=(atan2(UMEAN(I),VMEAN(I))+3.1415926)*57.29578
        if(DIRECTION(I).gt.360.)DIRECTION(I)=DIRECTION(I)-360.
        if(VSPEED.lt.5.)DIRECTION(I)=9999.
     enddo

     do K=2,MKX-1
        do I = ISTART,IEND
           DQ=.5*(Q(I,K+1)-Q(I,K-1))
           ! MCONV(I)=MCONV(I)+1.E5*OMEG(I,K)*DQ/G
           ! CONVERGENCIA DE UMIDADE DA COLUNA
           MCONV(I)=MCONV(I)+OMEG(I,K)*DQ/G 
           ! CONVERGENCIA DE UMIDADE DA COLUNA (OMEGA EM PA/S)
        enddo
     enddo
     do I = ISTART,IEND
       if(MCONV(I).lt. 0.)  MCONV(I) = 0.
     enddo

     !---  CUMULUS PARAMETERIZATION
     !srf- aqui se deve colocal o loop no ensemble dependente do tipo
     !     de cumulus
     !iens =1

     call CUP_enss(ngrid, mynum, m1, m2, m3, i0, j0, ipr, jpr,    &
          mgmxp, mgmyp, mgmzp, maxiens, maxens, maxens2,          &
          maxens3, ensdim, icoic, j, iens, ISTART, IEND,          &
          mix, mjx, mkx, massfln, massflx, iact_gr, iact_old_gr,  &
          xland, TER11, AA0, T, Q, TN, QO, PO, PRET, P, OUTT,     &
          OUTQ, OUTQC, DTIME, PSUR, US_Grell, VS_Grell, KDET, TCRIT, time,    &
          MCONV, OMEG, DIRECTION,                                 &
          !Insercao de variaveis para salvar  propriedades dos cumulus
          ierr4d(1,j,iens,ngrid), jmin4d(1,j,iens,ngrid),        &
          kdet4d(1,j,iens,ngrid), k224d(1,j,iens,ngrid),         &
          kbcon4d(1,j,iens,ngrid), ktop4d(1,j,iens,ngrid),       &
          kpbl4d(1,j,iens,ngrid),                                &
          kstabi4d(1,j,iens,ngrid), kstabm4d(1,j,iens,ngrid),    &
          xmb4d(1,j,iens,ngrid), edt4d(1,j,iens,ngrid),          &
          zcup5d(1,1,j,iens,ngrid),pcup5d(1,1,j,iens,ngrid),     &
          enup5d(1,1,j,iens,ngrid), endn5d(1,1,j,iens,ngrid),    &
          deup5d(1,1,j,iens,ngrid), dedn5d(1,1,j,iens,ngrid),    &
          zup5d(1,1,j,iens,ngrid), zdn5d(1,1,j,iens,ngrid),      &
          prup5d(1,1,j,iens,ngrid),clwup5d(1,1,j,iens,ngrid),    &
          tup5d(1,1,j,iens,ngrid),                               &
          !Insercao de variaveis para salvar  nas analises e  Cup_DIRECTION2
          upmf, dnmf, xierr, xktop, xkbcon, xk22, xjmin, xkdt,         &
          xiact_p, xiact_c)

     !srf-----print-------
     do I=ISTART,IEND
        if(j.eq.jpr   .and. i.eq.ipr) then

           do K=1,ktop4d(i,j,1,1)
              if(k.eq.1) then
                 write(2,'(a1,78a1)') ' ',(' ',m=1,78)

                 write(2,1113) ierr4d(i,j,iens,ngrid), time/3600.,  &
                      UPMF(I,J), DNMF(I,J), PRET(I)*3600.
                 write(2,1114) int(XIERR(i,J)), int(XKBCON(I,J)),   &
                      int(XKTOP(I,J)), int(XJMIN(I,J))
                 write(2,'(a1,78a1)') ' ',(' ',m=1,78)

              endif

              write(2,1115) K, outt(i,k)*86400., OUTQ(i,k)*86400.*1000.,  &
                   OUTQC(i,k)*86400.*1000.
1113          format(1x,i4,4E12.5)
1114          format(1x,4I4)
1115          format(1x,I4,3E12.5)

           enddo
        endif
     enddo

     !	  DO I=ISTART,IEND
     !	  if(ierr4d(i,j,iens,ngrid) .eq. 0) then
     !	  print*,'&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&'
     !	  print*,'2 mynum i j ierr=',mynum,i,j,ierr4d(i,j,iens,ngrid)
     !	  print*,'XMB=',xmb4d(i,j,iens,ngrid)
     !	  do  k=1,mkx
     !	  print*,k,zcup5d(k,i,j,iens,ngrid),deup5d(k,i,j,iens,ngrid)
     !	  ENDDO
     !	  print*,'&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&'
     !!          stop
     !	  endif
     !	  enddo
     !	  IF(J.EQ.J2) STOP
     !PRINT------------
     !

     !--- Output

     do K=1,MKX-1 
        kr = K + 1
        do I = ISTART,IEND
           !      Converte tendencia da temperatura (OUTT) em 
           !      tendencia de theta (OUTTEM)
           !      cp*T=Pi*Theta => cp dT/dt = Theta*dPi/dt + Pi*dTheta/dt,
           !      assumindo dPi/dt (=pt(kr,i,j)) << (exner/theta)*dTheta/dt:
           !      Exner's function = pp(kr,i,j)+pi0(kr,i,j)
           exner          = pp(kr,i,j) + pi0(kr,i,j)
           outtem(kr,i,j) = CP/exner   * OUTT(I,K) 
           ! tendencia do Theta  devida aos cumulus
           outrt(kr,i,j)  = OUTQ(I,K)  + OUTQC(I,K) 
           ! tendencia do Rtotal devida aos cumulus
        enddo
     enddo

     !srf-----print-------
     !      DO I = ISTART,IEND
     !         if(int(xierr(i,j)).eq.0) then
     !          write(mynum+4,122) mynum,i0+i,j0+j,xierr(i,j),pret(I)
     ! 122      format(1x,3i5,f13.4,f13.4)
     !       
     !       DO K=MKX-1,1,-1
     !       kr = K + 1
     !       write(mynum+4,121) mynum,i0+i,j0+j,kr,outt(i,kr)*86400.,pret(I)
     ! 121      format(1x,4i5,e13.4,f13.4)
     !          print*,'pr# i0+i j0+j prec=',mynum,i0+i,j0+j,PRET(I)
     !       ENDDO
     !         endif
     !      ENDDO
     !srf-----print-------

     do I = ISTART,IEND
        PRECIP(I,J)=PRET(I)

        !srf-----print-------
        !         if(PRECIP(I,J).GT.0.) then
        !         if(j.eq.2) then
        !          write(mynum+2,101) mynum,i0+i,j0+j,PRECIP(i,j)
        ! 101      format(1x,3i5,f13.4)
        !          print*,'pr# i0+i j0+j prec=',mynum,i0+i,j0+j,PRET(I)
        !         endif
        !srf-----print-------

     enddo

     do I = ISTART,IEND
        if(PRECIP(I,J).le.0.)then
           IACT_GR(I,J) =0
           PRECIP(I,J)  =0.
           do K=1,MKX
              kr = K + 1
              outtem(kr,i,j) = 0.
              ! comente para testes com conservacao (c0=0.)
              outrt(kr,i,j) =0.
              ! comente para testes com conservacao (c0=0.)
           enddo
           do K=1,ENSDIM
              MASSFLN(I,J,K)=0.
           enddo
        else
           IACT_GR(I,J)=1
        endif
     enddo

     !--- Salva nas analises

     do I = ISTART,IEND
        massflx(i,j) = dnmf(i,j)
        xiact_c(i,j) = float(IACT_GR(I,J))
        xiact_p(i,j) = float(iact_old_gr(I,J))
     enddo

  enddo     ! loop externo - j -

  !--- shallow convection

  return
end subroutine CUPARTH

!--------------------------------------------------------------------
!
!

subroutine CUP_enss(ngrid, mynum, m1, m2, m3, i0, j0, ipr, jpr,       &
     mgmxp, mgmyp, mgmzp, maxiens, maxens, maxens2, maxens3, ensdim,  &
     icoic, j, iens, ISTART, IEND, mix, mjx, mkx, massfln, massflx,   &
     iact_gr, iact_old_gr, xland, Z1, AAEQ, T, Q, TN, QO, PO, PRE, P, &
     OUTT, OUTQ, OUTQC, DTIME, PSUR, US, VS, KDET, TCRIT, time,       &
     mconv, omeg, direction, ierr4d, jmin4d, kdet4d, k224d, kbcon4d,  &
     ktop4d, kpbl4d, kstabi4d, kstabm4d, xmb4d, edt4d, zcup5d,        &
     pcup5d, enup5d, endn5d, deup5d, dedn5d, zup5d, zdn5d,            &
     !Lufla     p_lw5d,    &
     prup5d,clwup5d,tup5d, & !Lufla    &
     upmf, dnmf, xierr, xktop, xkbcon, xk22, xjmin, xkdt, xiact_p, xiact_c)

  ! USE Modules for Grell Parameterization
  use mem_scratch3_grell

  implicit none
  integer maxiens,maxens,maxens2,maxens3,ensdim
  integer mix,mjx,mkx,mgmxp, mgmyp, mgmzp
  integer nall,nens,iens,iedt,ktau
  integer nens3
  integer izero
  integer icoic
  integer, save :: ialloc
  data ialloc/0/
  real time

  !--- Input variables -----------------------------

  real mconv(mgmxp), Z1(mgmxp), direction(mgmxp), AAEQ(mgmxp),  &
       pre(mgmxp), PSUR(mgmxp)
  real T(mgmxp,mgmzp), Q(mgmxp,mgmzp), TN(mgmxp,mgmzp), QO(mgmxp,mgmzp),   &
       P(mgmxp,mgmzp), PO(mgmxp,mgmzp), US(mgmxp,mgmzp), VS(mgmxp,mgmzp),  &
       omeg(mgmxp,mgmzp)

  real massfln(mgmxp,mgmyp,ensdim)
  real massflx(mgmxp,mgmyp), xland(mgmxp,mgmyp)
  integer iact_gr(mgmxp,mgmyp), iact_old_gr(mgmxp,mgmyp), kdet(mgmxp)
  integer ngrid, mynum, i0, j0, ipr, jpr, m1, m2, m3

  !-------Salva parametros da CUP para uso no transporte convectivo:
  real  edt_average

  integer, dimension(mgmxp) :: ierr4d, jmin4d, kdet4d, k224d, kbcon4d,  &
       ktop4d, kpbl4d, kstabi4d, kstabm4d

  real, dimension(mgmxp) :: xmb4d, edt4d

  real, dimension(mgmzp,mgmxp) ::  &
       zcup5d,pcup5d, & !Lufla
       enup5d, endn5d, deup5d,  &
       dedn5d, zup5d, zdn5d, & !, p_lw5d Lufla
       prup5d,clwup5d,tup5d !Lufla

  !------Variables saved in RAMS Analisys
  ! use (m1,m2,m3) para dimensionar os vetores que sao
  ! escritos nas analises do RAMS
  real, dimension(m2,m3) :: upmf, dnmf, xierr, xktop, xkbcon, xjmin,    &
       xkdt, xiact_p, xiact_c,xk22

  !
  !--- Work variables - Allocatable in this point --
  !
  integer fquasi, fstab, fmconv, iresult
  integer ki, ip1, jp1, m
  integer I, J, K, ISTART, IEND, KK

  real mbdt
  !srf - added on 03-mar-2002

  !--- Output variables ----------------------------

  real OUTQC(mgmxp,mgmzp), OUTT(mgmxp,mgmzp), OUTQ(mgmxp,mgmzp)

  !--- Variables on cloud levels ------------------------------------

  real day, dz, tcrit, xl, pbcdif, outtes, dtime, outteq
  !srf_tmp
  real dellaqsum, dellaqcsum, dp
  !srf_tmp

  !--- New entrainment/detrainment related stuff --------------------

  real mentr_rate, mentrd_rate, entr_rate, radius, entrd_rate,  &
       massfld, zcutdown, edtmax, edtmin, depth_min, zkbmax,    &
       z_detr, zktop, dh, cap_maxs

  !----------------------End of memory allocation ----

  !
  !     if(ktau.gt.3.and.ktau.lt.7) ...
  !
  day=86400.
  !
  !--- specify entrainmentrate and detrainmentrate
  !
  !     radius=14000.-float(iens)*2000.
  radius=12000.
  !      radius=5000.
  fquasi=1
  fstab=0  
  fmconv=0
  !
  !--- gross entrainment rate (these may be changed later on in the
  !--- program, depending what your detrainment is!!)
  !
  entr_rate=.2/radius
  !
  !--- entrainment of mass
  !
  !      mentrd_rate=0.
  mentrd_rate=entr_rate
  mentr_rate =entr_rate
  !
  !--- initial detrainmentrates
  !
  do k=1,mkx
     do i=istart,iend
        cd(i,k)  = 0.1*entr_rate
        !cd(i,k) = (0.5-float(iens)*.1)*entr_rate
        !cd(i,k) = 0.
        cdd(i,k) = 0.
     enddo
  enddo
  !
  !--- max/min allowed value for epsilon
  !    (ratio downdraft base mass flux/updraft
  !    base mass flux
  !
  edtmax=.95
  edtmin=.2
  !
  !--- minimum depth (m), clouds must have
  !
  depth_min=500.
  !
  !--- maximum depth (mb) of capping 
  !--- inversion (larger cap = no convection)
  !
  !     cap_maxs=75.
  !     if(iens.eq.2)then
  !        cap_maxs=200.
  !     endif
  !     if(iens.eq.3)then
  !          radius=8.
  !     endif
  !     if(iens.eq.4)then
  !          do k=1,mkx
  !          do i=istart,iend
  !            cd(i,k)=0.
  !            cd(i,k)=0.1*entr_rate
  !          enddo
  !          enddo
  !     endif
  !
  !---- for kbcon_cin test
  cap_maxs=10.
  if (iens.eq.2) cap_maxs=20.
  if (iens.eq.3) cap_maxs=30.
  if (iens.eq.4) cap_maxs=40.
!  cap_maxs = 20.
  cap_maxs = 100.

  do I=ISTART,IEND
     !  DO 7 I=ISTART,IEND
     aa0(i)=0.
     aa1(i)=0.
     aad(i)=0.
     kstabm(i)=mkx-2
     if (aaeq(i).lt.0.) then
        ierr(i)=20
     else
        IERR(i)=0
        XIERR(i,j)=0.
        cwf(i,j)=0.
        pwf(i,j)=0.
        pwdf(i,j)=0.
        eddt(i,j)=0.
        !         xktop(i,j)=0.
        !         xkbas(i,j)=0.
        xmass(i,j)=0.
        predb(i,j)=0.
     endif
     !7    CONTINUE
  enddo

  !--- first check for upstream convection

  do i=istart,iend
     if (ierr(i).eq.0) then
        iresult=0
        massfld=0.
        !srf  CALL cup_direction2(i,j,direction,iact_old_gr,mix,mjx,  &
        !srf  mgmxp,mgmyp,massflx,iresult,ensdim,0,0,maxens3,massfld)

        cap_max(i)=cap_maxs
        if (iresult.eq.1) then
           cap_max(i)=cap_maxs+20.
        endif
     endif
  enddo

  !--- max height(m) above ground where updraft air can originate

  zkbmax=4000.

  !--- height(m) above which no downdrafts are allowed to originate

  zcutdown=3000.

  !--- depth(m) over which downdraft detrains all its mass

  z_detr=1250.

  do nens=1,maxens
     mbdt_ens(nens)=(float(nens)-3.)*dtime*1.e-3+dtime*5.E-03
     !        mbdt_ens(nens)=(float(nens)-1.5)*dtime*2.e-3+dtime*5.E-03
  enddo
  do nens=1,maxens2
     !        edt_ens(nens)=.7-float(nens)*.1
     edt_ens(nens)=.95-float(nens)*.01
  enddo
  !     if(j.eq.jpr)then
  !       print *,'radius ensemble ',iens,radius
  !       print *,mbdt_ens
  !       print *,edt_ens
  !     endif
  !
  !--- environmental conditions, FIRST HEIGHTS
  !

  do i=istart,iend
     if(ierr(i).ne.20)then
        do k=1,maxens*maxens2*maxens3
           xf_ens(  i,j,(iens-1)*maxens*maxens2*maxens3+k)= 0.
           pr_ens(  i,j,(iens-1)*maxens*maxens2*maxens3+k)= 0.
           outt_ens(i,j,(iens-1)*maxens*maxens2*maxens3+k)= 0.
        enddo
     endif
  enddo

  !--- calculate moist static energy, heights, qes
  !
  call cup_env(j, ipr, jpr, z, qes, he, hes, t, q, p, z1, mix, mgmxp,  &
       mkx, mgmzp, istart, iend, psur, ierr, tcrit, 0)
  call cup_env(j, ipr, jpr, zo, qeso, heo, heso, tn, qo, po, z1, mix,  &
       mgmxp, mkx, mgmzp, istart, iend, psur, ierr, tcrit, 0)
  !
  !--- environmental values on cloud levels
  !
  call cup_env_clev(j, ipr, jpr, t, qes, q, he, hes, z, p, qes_cup,  &
       q_cup, he_cup, hes_cup, z_cup, p_cup, gamma_cup, t_cup, psur, &
       mix, mgmxp, mkx, mgmzp, istart, iend, ierr, z1)
  call cup_env_clev(j, ipr, jpr, tn, qeso, qo, heo, heso, zo, po,    &
       qeso_cup, qo_cup, heo_cup, heso_cup, zo_cup, po_cup,          &
       gammao_cup, tn_cup, psur, mix, mgmxp, mkx, mgmzp, istart,     &
       iend, ierr, z1)

  !srf
  !       if(j.eq.jpr) then
  !        i=ipr
  !        do k=1,mkx
  !        write(6,'(3i4,5F12.4)')i,j,k,t(i,k),tn(i,k),he(i,k),
  !     &                          he_cup(i,k)
  !        enddo
  !        stop
  !       endif
  !srf

  do i=istart,iend
     if (ierr(i).eq.0) then

        do k=1,mkx
           if (zo_cup(i,k).gt.zkbmax+z1(i)) then
              kbmax(i)=k
              !GO TO 25
              exit
           endif
        enddo
!25      CONTINUE

        !--- level where detrainment for downdraft starts

        do k=1,mkx
           if (zo_cup(i,k).gt.z_detr+z1(i)) then
              kdet(i)=k
              !GO TO 26
              exit
           endif
        enddo
!26      CONTINUE

     endif
  enddo

  !--- DETERMINE LEVEL WITH HIGHEST MOIST STATIC ENERGY CONTENT - K22

  call MAXIMI(HEO_CUP, mix, mgmxp, mkx, mgmzp, 3, KBMAX, K22,  &
       ISTART, IEND, ierr)
  do I=ISTART,IEND
     !  DO 36 I=ISTART,IEND
     if (ierr(I).eq.0.) then
        if (K22(I).ge.KBMAX(i)) ierr(i)=2
     endif
     !36   CONTINUE
  enddo

  !--- DETERMINE THE LEVEL OF CONVECTIVE CLOUD BASE  - KBCON

  !Grell sugested to test cup_kbcon -18dec2001

  call cup_kbcon(1, k22, kbcon, heo_cup, heso_cup, mix, mgmxp,  &
       mkx, mgmzp, istart, iend, ierr, kbmax, po_cup, cap_max)

  !srf 30-jan-2002. This routine blow up the model
  !CALL cup_kbcon_cin(1, k22,kbcon, heo_cup, heso_cup, z, tn_cup,  &
  !     qeso_cup, mix, mgmxp, mkx, mgmzp, istart, iend, ierr,      &
  !     kbmax, po_cup, cap_max)

  !srf 
  !      DO I=ISTART,IEND
  !       if(i.eq.ipr.and.j.eq.jpr) then
  !       print*,'-------------------------------'
  !       print*,'MYNUM I J=',mynum,i,j
  !       print*,'k22-kbcon-ierr=',k22(i),kbcon(i),ierr(i)
  !       print*,'-------------------------------'
  !       endif
  !      enddo
  !srf

  !--- Increase detrainment in stable layers

  call MINIMI(HEso_cup, mix, mgmxp, mkx, mgmzp, Kbcon, kstabm,  &
       kstabi, ISTART, IEND, ierr)
  !Modificated by Grell - 07dec2001
  !      DO I=ISTART,IEND
  !      IF(ierr(I).eq.0.)THEN
  !        if(kstabm(i)-1.gt.kstabi(i))then
  !           do k=kstabi(i),kstabm(i)-1
  !            cd(i,k)=cd(i,k-1)+0.5*entr_rate
  !     1               /float(kstabm(i)-kstabi(i))
  !             if(cd(i,k).gt.0.7*entr_rate)cd(i,k)=0.7*entr_rate
  !           enddo
  !        ENDIF
  !      ENDIF
  !      ENDDO
  do I=ISTART,IEND
     if (ierr(I).eq.0.) then
        if (kstabm(i)-1.gt.kstabi(i)) then
           do k=kstabi(i),kstabm(i)-1
              cd(i,k)=cd(i,k-1)+1.5*entr_rate
              if (iens.gt.4) then
                 cd(i,k) = cd(i,k-1)+float(iens-4)*entr_rate/  &
                      float(kstabm(i)-kstabi(i))
              else
                 cd(i,k)=cd(i,k)
              endif
              if (cd(i,k).gt.10.0*entr_rate) cd(i,k)=10.0*entr_rate
           enddo
        endif
     endif
  enddo

  !--- Calculate incloud moist static energy

  call cup_up_he(k22, hkb, z_cup, cd, mentr_rate, he_cup, hc,   &
       mix, mgmxp, mkx, mgmzp, kbcon, ierr, istart, iend, dby,  &
       he, hes_cup)
  call cup_up_he(k22, hkbo, zo_cup, cd, mentr_rate, heo_cup, hco,  &
       mix, mgmxp, mkx, mgmzp, kbcon, ierr, istart, iend, dbyo,    &
       heo, heso_cup)

  !--- Determine cloud top - KTOP
  !
  !srf--print--------------------
  do I=ISTART,IEND
     if(i.eq.ipr.and.j.eq.jpr) then
        print*,'-----------K T O P-1-----------'
        print*,'MYNUM I J=',mynum,i,j
        print*,'k22-kbcon-ierr-Ktop=',k22(i),kbcon(i),ierr(i),ktop(i)
        do K=KBCON(I)+1,MKX-2
           print*,'k-heo_cup-dbyo=',k,heo_cup(i,k),dbyo(i,k)
        enddo
        print*,'-------------------------------'
     endif
  enddo
  !srf--print--------------------

  !srf
  !IF(j.EQ.jpr)  &
  !     CALL cup_ktop(1, dbyo, kbcon, ktop, mix, mgmxp, mkx, mgmzp,  &
  !     istart, iend, ierr)
  !srf
  call cup_ktop(1, dbyo, kbcon, ktop, mix, mgmxp, mkx, mgmzp, istart, &
       iend, ierr)

  !srf--print--------------------
  !DO I=ISTART,IEND
  !   if(i.eq.ipr.and.j.eq.jpr) then
  !     print*,'-----------K T O P--2----------'
  !     print*,'MYNUM I J=',mynum,i,j
  !     print*,'k22-kbcon-ierr-Ktop=',k22(i),kbcon(i),ierr(i),ktop(i)
  !     do K=KBCON(I)+1,MKX-2
  !        print*,'k-heo_cup-dbyo=',k,heo_cup(i,k),dbyo(i,k)
  !     enddo
  !     print*,'-------------------------------'
  !   endif
  !enddo
  !srf--print--------------------

  do I=ISTART,IEND
     ! DO 37 I=ISTART,IEND
     kzdown(i)=0
     if (ierr(i).eq.0) then
        zktop = (zo_cup(i,ktop(i))-z1(i))*.6
        zktop = min(zktop+z1(i),zcutdown+z1(i))
        do k=1,mkx
           if (zo_cup(i,k).gt.zktop) then
              kzdown(i)=k
              !GO TO 37
              exit
           endif
        enddo
     endif
!37   CONTINUE
  enddo

  !--- Downdraft originating level - JMIN

  call MINIMI(HEso_cup, mix, mgmxp, mkx, mgmzp, K22, kzdown,  &
       JMIN, ISTART, IEND, ierr)

  !srf--print-------------------
  !      DO I=ISTART,IEND
  !       if(i.eq.ipr.and.j.eq.jpr) then
  !       print*,'-------------------------------'
  !          write(mynum+2,111) mynum,i0+i,j0+j,Jmin(i),Ktop(i)
  ! 111      format(1x,5i5)
  !       print*,'MYNUM I J=',mynum,i,j
  !       print*,'k22-kbcon-ierr=',k22(i),kbcon(i),ierr(i)
  !       print*,'ktop-he_cup-dbyo=',ktop(i),he_cup(i,12),dbyo(i,12)
  !       print*,'Ktop Jmin=',ktop(i),Jmin(i)
  !       print*,'-------------------------------'
  !       endif
  !      enddo
  !srf--print--------------------

  do I=ISTART,IEND
     !  DO 100 I=ISTART,IEND
     if (ierr(I).eq.0.) then

        !--- Check whether it would have buoyancy, if there where
        !--- no entrainment/detrainment

101     continue
        if (jmin(i)-1.lt.KDET(I)) kdet(i)=jmin(i)-1
        if (jmin(i).ge.Ktop(I)-1) jmin(i)=ktop(i)-2
        ki=jmin(i)

        !srf--print-------------------
        !       if(i.eq.ipr.and.j.eq.jpr) then
        !        write(mynum+2,111) mynum,i0+i,j0+j,Jmin(i),Ktop(i)
        !        print*,'-------------------------------'
        !        print*,'MYNUM I J=',mynum,i,j
        !        print*,'Ktop Jmin=',ktop(i),Jmin(i)
        !        print*,'-------------------------------'
        !       endif
        !srf--print--------------------

        !     hcdo(i,ki)=heo_cup(i,ki)
        hcdo(i,ki) = heso_cup(i,ki)
        DZ         = Zo_cup(i,Ki+1)-Zo_cup(i,Ki)
        dh         = dz*(HCDo(i,Ki)-heso_cup(i,ki))
        dh         = 0.

        do k=ki-1,1,-1
           !        hcdo(i,k)=heo_cup(i,jmin(i))
           hcdo(i,k) = heso_cup(i,jmin(i))
           DZ        = Zo_cup(i,K+1)-Zo_cup(i,K)
           dh        = dh+dz*(HCDo(i,K)-heso_cup(i,k))
           if (dh.gt.0.) then
              jmin(i)=jmin(i)-1
              if (jmin(i).gt.3) then
                 GO TO 101
              else if (jmin(i).le.3) then
                 ierr(i)=9
                 !CYCLE  !GO TO 100
                 GO TO 100
              endif
           endif
        enddo

        if (JMIN(I).le.3) then
           ierr(i)=4
        endif

     endif
100  continue
  enddo

  !srf 
  !      DO I=ISTART,IEND
  !       if(i.eq.ipr.and.j.eq.jpr)
  !     &  print*,'ktop-JMIN-ierr=',ktop(i),jmin(i),ierr(i)
  !      enddo
  !srf

  !--- Must have at least depth_min m between cloud convective
  !    base and cloud top.

  do i=istart,iend
     if (ierr(I).eq.0.) then
        if (-zo_cup(I,KBCON(I))+zo_cup(I,KTOP(I)).lt.  &
             depth_min) then
           ierr(i)=6
        endif
     endif
  enddo

  !--- Normalized updraft mass flux profile

  call cup_up_nms(zu, z_cup, mentr_rate, cd, kbcon, ktop,  &
       mix, mgmxp, mkx, mgmzp, istart, iend, ierr, k22)
  call cup_up_nms(zuo, zo_cup, mentr_rate, cd, kbcon, ktop,  &
       mix, mgmxp, mkx, mgmzp, istart, iend, ierr, k22)

  !--- Normalized downdraft mass flux profile,also work on
  !    bottom detrainment
  !--- in this routine

  call cup_dd_nms(zd, z_cup, cdd, mentrd_rate, jmin, ierr,    &
       mix, mgmxp, mkx, mgmzp, istart, iend, 0, kdet, z1)
  call cup_dd_nms(zdo, zo_cup, cdd, mentrd_rate, jmin, ierr,  &
       mix, mgmxp, mkx, mgmzp, istart, iend, 1, kdet, z1)

  !--- Downdraft moist static energy

  call cup_dd_he(hes_cup, zd, hcd, z_cup, cdd, mentrd_rate,  &
       jmin, ierr, mix, mgmxp, mkx, mgmzp, istart, iend,he,  &
       kdet, dbyd, he_cup)
  call cup_dd_he(heso_cup, zdo, hcdo, zo_cup, cdd,           &
       mentrd_rate, jmin, ierr, mix, mgmxp, mkx, mgmzp,      &
       istart, iend,heo, kdet, dbydo, he_cup)

  !--- Calculate moisture properties of downdraft

  call cup_dd_moisture(j, zd, hcd, hes_cup, qcd, qes_cup,  &
       pwd, q_cup, z_cup, cdd, mentrd_rate, jmin, ierr,    &
       gamma_cup, pwev, mix, mgmxp, mkx, mgmzp, istart,    &
       iend, bu, qrcd, q, he, hc, t_cup, 2)
  call cup_dd_moisture(j, zdo, hcdo, heso_cup, qcdo,       &
       qeso_cup, pwdo, qo_cup, zo_cup, cdd, mentrd_rate,   &
       jmin, ierr, gammao_cup, pwevo, mix, mgmxp, mkx,     &
       mgmzp, istart, iend, bu, qrcdo, qo, heo, hco, tn_cup, 1)

  !--- Calculate moisture properties of updraft

  call cup_up_moisture(ierr, z_cup, qc, qrc, pw, pwav,     &
       kbcon, ktop, mix, mgmxp, mkx, mgmzp, istart, iend,  &
       cd, dby, mentr_rate, q, GAMMA_cup, zu, qes_cup,     &
       k22, q_cup)
  call cup_up_moisture(ierr, zo_cup, qco, qrco, pwo, pwavo,  &
       kbcon, ktop, mix, mgmxp, mkx, mgmzp, istart, iend,    &
       cd, dbyo,                                             &
       mentr_rate, qo, GAMMAo_cup, zuo, qeso_cup, k22, qo_cup)

  !--- Calculate workfunctions for updrafts

  call cup_up_aa0(aa0, z, zu, dby, GAMMA_CUP, t_cup, kbcon,  &
       ktop, mix, mgmxp, mkx, mgmzp, istart, iend, ierr)
  call cup_up_aa0(aa1, zo, zuo, dbyo, GAMMAo_CUP, tn_cup,    &
       kbcon, ktop, mix, mgmxp, mkx, mgmzp, istart, iend, ierr)
  do i=istart,iend
     if (ierr(i).eq.0) then
        if (aa1(i).eq.0.) then
           ierr(i)=17
        endif
     endif
  enddo

  !--- Determine downdraft strength in terms of windshear

  ! Init. array EDTC with Zeros - ALF

  do k=1,maxens2
     do i=1,mgmxp
        EDTC(i,k) = 0.
     enddo
  enddo

  call cup_dd_edt(ierr, us, vs, zo, ktop, kbcon, edt, po,   &
       pwavo, pwevo, mix, mgmxp, mkx, mgmzp, istart, iend,  &
       edtmax, edtmin, maxens2, edtc, vshear, sdp, vws)

  !srf - Big loop starts here!

  do iedt=1,maxens2
     !DO 250 iedt=1,maxens2
     do i=istart,iend
        if (ierr(i).eq.0) then
           edt(i)  = edtc(i,iedt)
           edto(i) = edtc(i,iedt)
           edtx(i) = edtc(i,iedt)
        endif
     enddo
     do k=1,mkx
        do i=istart,iend
           dellat_ens(i,k,iedt) = 0.
           dellaq_ens(i,k,iedt) = 0.
           dellaqc_ens(i,k,iedt)= 0.
           pwo_ens(i,k,iedt)    = 0.
        enddo
     enddo

     if (j.eq.jpr.and.iedt.eq.1000000) then
        i=ipr
        print *,'in 250 loop ', iedt, edt(ipr)
        if (ierr(i).eq.0.or.ierr(i).eq.3) then
           print *, k22(I), kbcon(i), ktop(i), jmin(i)
           print *, edt(i), aa0(i), aa1(i)
           do k=1,mkx
              print *,z(i,k),he(i,k),hes(i,k)
           enddo
           do k=1,ktop(i)+1
              print *,zu(i,k),zd(i,k),pw(i,k),pwd(i,k)
           enddo
        endif
     endif
     do I=ISTART,IEND
        aad(i)=0.
     enddo
     !     DO I=ISTART,IEND
     !       if(ierr(i).eq.0)then
     !        eddt(i,j)=edt(i)
     !        EDTX(I)=EDT(I)
     !        BU(I)=0.
     !        BUO(I)=0.
     !       endif
     !     enddo
     !
     !--- Downdraft workfunctions

     !CALL cup_dd_aa0(edt, ierr, aa0, jmin, gamma_cup, t_cup,  &
     !     hcd, hes_cup, z, mix, mgmxp, mkx, mgmzp, istart,    &
     !     iend, zd)

     call cup_dd_aa0(edto, ierr, aad, jmin, gammao_cup, tn_cup, &
          hcdo, heso_cup, zo, mix, mgmxp, mkx, mgmzp, istart,   &
          iend, zdo)

     !--- Change per unit mass that a model cloud would
     !    modify the environment

     !--- 1. in bottom layer

     !srf trocando 'po'  por   'po_cup'
     !call cup_dellabot(ipr,jpr,heo_cup,ierr,zo_cup,po,hcdo,edto,
     !zdo,cdd,heo,mix,mgmxp,mkx,mgmzp,istart,iend,dellah,1,j,mentrd_rate,zo)
     !call cup_dellabot(ipr,jpr,qo_cup,ierr,zo_cup,po,qrcdo,edto,
     !zdo,cdd,qo,mix,mgmxp,mkx,mgmzp,istart,iend,dellaq,2,j,mentrd_rate,zo)

     call cup_dellabot(0, 0, heo_cup, ierr, zo_cup, po_cup,   &
          hcdo, edto, zdo, cdd, heo, mix, mgmxp, mkx, mgmzp,  &
          istart, iend, dellah, 1, j, mentrd_rate, zo)
     call cup_dellabot(ipr, jpr, qo_cup, ierr, zo_cup,    &
          po_cup, qrcdo, edto, zdo, cdd, qo, mix, mgmxp,  &
          mkx, mgmzp, istart, iend, dellaq, 2, j,         &
          mentrd_rate, zo)
     !srf trocando 'po'  por   'po_cup'

     !--- 2. everywhere else

     call cup_dellas(ierr, zo_cup, po_cup, hcdo, edto, zdo,  &
          cdd, heo, mix, mgmxp, mkx, mgmzp, istart, iend,    &
          dellah, 1, j, mentrd_rate, zuo, cd, hco, ktop,     &
          k22, kbcon, mentr_rate, jmin, heo_cup, kdet, k22,  &
          0, 0, 'deep')

     !-- Take out cloud liquid water for detrainment

     do k=1,mkx
        do i=istart,iend
           scr1(i,k)=0.
           dellaqc(i,k)=0.
           if (ierr(i).eq.0) then
              scr1(i,k)=qco(i,k)-qrco(i,k)

              if (k.eq.ktop(i)-0)                    &
                   dellaqc(i,k)=.01*zuo(i,ktop(i))*  &
                   qrco(i,ktop(i))*9.81/(po_cup(i,k)-po_cup(i,k+1))

              if (k.lt.ktop(i).and.k.gt.kbcon(i)) then
                 dz = zo_cup(i,k+1)-zo_cup(i,k)
                 dellaqc(i,k) = .01*9.81*cd(i,k)*dz*zuo(i,k)*  &
                      .5*(qrco(i,k)+qrco(i,k+1))/              &
                      (po_cup(i,k  )-po_cup(i,k+1))
              endif
           endif
        enddo
     enddo

     call cup_dellas(ierr, zo_cup, po_cup, qrcdo, edto, zdo,  &
          cdd, qo, mix, mgmxp, mkx, mgmzp, istart, iend,      &
          dellaq, 2, j, mentrd_rate, zuo, cd, scr1, ktop,     &
          k22, kbcon, mentr_rate, jmin, qo_cup, kdet, k22,    &
          ipr, jpr, 'deep')

     !srf-----print-------
     do i=istart,iend
        if (ierr(i).eq.0) then
           if (j.eq.jpr.and.i.eq.ipr) then
              dellaqsum  =0.
              dellaqcsum =0.
              do k=1,mkx
                 dp         = -100.*(p_cup(i,k+1)-p_cup(i,k))
                 !	   dp=100.*(po(i,k-1)-po(i,k))
                 !	   if(k.eq.1) dp=996.2952
                 dellaqsum  = dellaqsum + dellaq(i,k) * dp/9.81
                 dellaqcsum = dellaqcsum + dellaqc(i,k)* dp/9.81

                 if (k.eq.1) then
                    write(6,'(a1,78a1)') ' ',('-',m=1,78)
                    print *, 'i k Z_cup P_cup dp dellaq dellaq',  &
                         ' dellaqsum dellaqcsum'
                 endif
                 write(6,'(2i4,5F12.4,3f16.4)') i, k,           &
                      Z_cup(i,k),p_cup(i,k), dp,                &
                      1.e+3*86400.*dellaq(i,k),                 &
                      1.e+3*86400.*dellaqc(i,k),                &
                      1.e+3*86400.*dellaqsum,                   &
                      1.e+3*86400.*dellaqcsum,                  &
                      1.e+3*86400.*(dellaqsum+dellaqcsum)
                 !& 100.*(p_cup(i,k-1)-p_cup(i,k))/
                 !& (g*(z_cup(i,k)-z_cup(i,k-1)))
                 ! densidade
              enddo
           endif
        endif
     enddo
     !srf-----print-------

     !--- Using dellas, calculate changed environmental profiles

     !Changed by Grell - 07dec2001
     !      do 200 nens=1,maxens
     !      mbdt=mbdt_ens(nens)
     !     do i=istart,iend 
     !      xaa0_ens(i,nens)=0.
     !      enddo
     !      do k=1,mkx-1
     !      do i=istart,iend
     !         if(ierr(i).eq.0)then
     !            XHE(I,K)=DELLAH(I,K)*MBDT+HEO(I,K)
     !            XQ(I,K)=DELLAQ(I,K)*MBDT+QO(I,K)
     !            DELLAT(I,K)=(1./1004.)*  &
     !                (DELLAH(I,K)-2.5E06*DELLAQ(I,K))
     !            XT_Grell(I,K)= DELLAT(I,K)*MBDT+TN(I,K)
     !            IF(XQ(I,K).LE.0.)XQ(I,K)=1.E-08
     !         ENDIF
     !      enddo
     !      enddo
     do nens=1,maxens
        !DO 200 nens=1,maxens
        mbdt=mbdt_ens(nens)
        do i=istart,iend
           xaa0_ens(i,nens)=0.
        enddo
        do k=1,mkx-1
           do i=istart,iend
              dellat(i,k)=0.
              if (ierr(i).eq.0) then
                 XHE(I,K)   = DELLAH(I,K)*MBDT + HEO(I,K)
                 XQ(I,K)    = DELLAQ(I,K)*MBDT +  QO(I,K)
                 DELLAT(I,K)= (1./1004.)*(DELLAH(I,K)-2.5E06*DELLAQ(I,K))
                 XT_Grell(I,K)    = DELLAT(I,K)*MBDT +  TN(I,K)
                 if (XQ(I,K).le.0.) XQ(I,K)=1.E-08
                 if (i.eq.ipr.and.j.eq.jpr) then
                    print *, k, DELLAH(I,K), DELLAQ(I,K), DELLAT(I,K)
                 endif
              endif
           enddo
        enddo
        !
        do i=istart,iend
           if (ierr(i).eq.0) then
              XHE(I,mkx) = HEO(I,mkx)
              XQ(I,mkx)  = QO(I,mkx)
              XT_Grell(I,mkx)  = TN(I,mkx)
              if (XQ(I,mkx).le.0.) XQ(I,mkx) = 1.E-08
           endif
        enddo

        !--- Calculate moist static energy, heights, qes

        call cup_env(j, ipr, jpr, xz, xqes, xhe, xhes, xt_grell,  &
             xq, po, z1, mix, mgmxp, mkx, mgmzp, istart,    &
             iend, psur, ierr, tcrit, 2)

        !--- Environmental values on cloud levels

        call cup_env_clev(j, ipr, jpr, xt_grell, xqes, xq, xhe,   &
             xhes, xz, po, xqes_cup, xq_cup, xhe_cup,       &
             xhes_cup, xz_cup, po_cup, gamma_cup, xt_cup,   &
             psur, mix, mgmxp, mkx, mgmzp, istart, iend,    &
             ierr, z1)

        !**************************** Static Control

        !--- Moist static energy inside cloud

        do i=istart,iend
           if (ierr(i).eq.0) then
              xhkb(i)=xhe(i,k22(i))
           endif
        enddo
        call cup_up_he(k22, xhkb, xz_cup, cd, mentr_rate,  &
             xhe_cup, xhc, mix, mgmxp, mkx, mgmzp, kbcon,  &
             ierr, istart, iend, xdby, xhe, xhes_cup)

        !--- Normalized mass flux profile

        call cup_up_nms(xzu, xz_cup, mentr_rate, cd, kbcon,  &
             ktop, mix, mgmxp, mkx, mgmzp, istart, iend,     &
             ierr, k22)
        call cup_dd_nms(xzd, xz_cup, cdd, mentrd_rate, jmin, &
             ierr, mix, mgmxp, mkx, mgmzp, istart, iend, 1,  &
             kdet, z1)

        !--- Moisture downdraft

        call cup_dd_he(xhes_cup, xzd, xhcd, xz_cup, cdd,     &
             mentrd_rate, jmin, ierr, mix, mgmxp, mkx,       &
             mgmzp, istart, iend, xhe, kdet, dbyd, xhe_cup)
        call cup_dd_moisture(j, xzd, xhcd, xhes_cup, xqcd,   &
             xqes_cup, xpwd, xq_cup, xz_cup, cdd,            &
             mentrd_rate, jmin, ierr, gamma_cup, xpwev,mix,  &
             mgmxp, mkx, mgmzp, istart, iend, bu, xqrcd, xq, &
             xhe, xhc, xt_cup,3)

        !--- Moisture updraft

        call cup_up_moisture(ierr, xz_cup, xqc, xqrc, xpw,   &
             xpwav, kbcon, ktop, mix, mgmxp, mkx, mgmzp,     &
             istart, iend, cd, xdby, mentr_rate, xq,         &
             GAMMA_cup, xzu, xqes_cup, k22, xq_cup)
        !
        !--- Workfunctions for updraft
        !
        call cup_up_aa0(xaa0, xz, xzu, xdby, GAMMA_CUP,   &
             xt_cup, kbcon, ktop, mix, mgmxp, mkx, mgmzp, &
             istart, iend,ierr)

        !--- Workfunctions for downdraft
        !
        !
        do i=istart,iend 
           if (ierr(i).eq.0) then
              xaa0_ens(i,nens) = xaa0(i)
              nall = (iens-1)*maxens3*maxens*maxens2 +  &
                   (iedt-1)*maxens*maxens3 + (nens-1)*maxens3

              do k=1,mkx
                 if (k.le.ktop(i)) then
                    do nens3=1,maxens3
                       if (nens3.eq.7) then
                          !--- b=0
                          pr_ens(i,j,nall+nens3) = pr_ens(i,j,nall+nens3) +  &
                               pwo(i,k)+edto(i)*pwdo(i,k)

                          !--- b=beta
                       else if (nens3.eq.8) then
                          pr_ens(i,j,nall+nens3) = pr_ens(i,j,nall+nens3) +  &
                               pwo(i,k)
                          !--- b=beta/2
                       else if (nens3.eq.9) then
                          pr_ens(i,j,nall+nens3) = pr_ens(i,j,nall+nens3) +  &
                               pwo(i,k)+.5*edto(i)*pwdo(i,k)
                       else
                          pr_ens(i,j,nall+nens3) =  pr_ens(i,j,nall+nens3) +  &
                               pwo(i,k)+edto(i)*pwdo(i,k)
                       endif
                    enddo
                 endif
              enddo

              do nens3=1,maxens3
                 outt_ens(i,j,nall+nens3) = dellat(i,1)
              enddo
           endif
        enddo
     !200     CONTINUE
     enddo

     !--- LARGE SCALE FORCING

     !------- CHECK wether aa0 should have been zero

     call MAXIMI(HE_CUP, mix, mgmxp, mkx, mgmzp, 3,   &
          KBMAX, K22x, ISTART, IEND, ierr)
     do I=ISTART,IEND
        if (ierr(i).eq.0) then
           if (K22x(I).ge.KBMAX(i)) ierr(i)=998
        endif
     enddo

     !--- DETERMINE THE LEVEL OF CONVECTIVE CLOUD BASE - KBCON

     !CALL cup_kbcon(2,k22x,kbconx,he_cup,hes_cup,
     !1    mix,mgmxp,mkx,mgmzp,istart,iend,ierr,kbmax,p_cup,cap_max)

!!!! Commented just for while 
    !Grell sugested on 07-dec-2001 to comment the 2 calls below in order
     !to start to test the ensemble:
     !call cup_kbcon_cin(1,k22x,kbconx,he_cup,hes_cup,z,t_cup,
     !1   qes_cup,mix,mgmxp,mkx,mgmzp,istart,iend,ierr,kbmax,p_cup,cap_max)
     !
     !call cup_ktop(2,dby,kbconx,ktopx,mix,mgmxp,mkx,mgmzp,istart,iend,ierr)
!!!!!
!srf - 31-jan-2003 - ensemble of closures
!     CALL cup_forcing_ens1(aa0, aa1, xaa0_ens, mbdt_ens,  &
!          dtime, xmb, ierr, mix, mgmxp, mjx, mgmyp, mkx,  &
!          mgmzp, istart, iend, xf_ens, j, fquasi, fstab,  &
!          'deeps', xland, maxens, iens, iedt, maxens2,    &
!          ipr, jpr, maxens3, mconv, omeg, zdo, kbcon,     &
!          zuo, pr_ens, edto, aad, kdet, massflx,          &
!          iact_old_gr, direction, ensdim, massfln,        &
!          xff_ens3, xk)

     call cup_forcing_ens_16(aa0,aa1,xaa0_ens,mbdt_ens,dtime,          &
          xmb,ierr,mix,mgmxp,mjx,mgmyp,mkx,mgmzp,istart,iend,xf_ens,j, &
          fquasi,fstab,'deeps',xland,maxens,iens,iedt,maxens2,ipr,jpr, &
          maxens3,mconv,omeg,zdo,kbcon,zuo,pr_ens,edto,aad,kbcon,       &
          massflx,iact_old_gr,direction,ensdim,                        &
          massfln,massfld,iresult,xff_ens3, xk,p_cup,ktop,icoic)       

     !
     !       DO I=ISTART,IEND
     !          xierr(i,j)=float(ierr(i))
     !          IF(ierr(i).eq.0)then
     !             cwf(i,j)=aa0(i)
     !             pwf(i,j)=pwav(i)
     !             pwdf(i,j)=pwev(i)
     !             xmass(i,j)=xmb(i)
     !           elseif(ierr(i).ne.0.and.ierr(i).ne.20)then
     !             cwf(i,j)=0.
     !             pwf(i,j)=0.
     !             pwdf(i,j)=0.
     !             xmass(i,j)=0.
     !          endif
     !       enddo
     !
     do k=1,mkx
        do i=istart,iend
           if (ierr(i).eq.0) then
              dellat_ens(i,k,iedt)  =  dellat(i,k)
              dellaq_ens(i,k,iedt)  =  dellaq(i,k)
              dellaqc_ens(i,k,iedt) = dellaqc(i,k)
              pwo_ens(i,k,iedt)  = pwo(i,k)+edt(i)*pwdo(i,k)
           else 
              dellat_ens(i,k,iedt)  = 0.
              dellaq_ens(i,k,iedt)  = 0.
              dellaqc_ens(i,k,iedt) = 0.
              pwo_ens(i,k,iedt)     = 0.
           endif
        enddo
     enddo
  !250  CONTINUE
  enddo

  !--- FEEDBACK

!  CALL cup_output_ens(xf_ens, ierr, dellat_ens, dellaq_ens, dellaqc_ens,   &
!       outt, outq, outqc, pre, pwo_ens, xmb, ktop, mix, mgmxp, mjx, mgmyp, &
!       mkx, mgmzp, istart, iend, j, 'deep', maxens2, maxens, ipr, jpr,     &
!       iens, pr_ens, outt_ens, maxens3, ensdim, massfln, xfac1)
              
  call cup_output_ens(xf_ens,ierr,dellat_ens,dellaq_ens,dellaqc_ens,   &
       outt,outq,outqc,pre,pwo_ens,xmb,ktop,mix,mgmxp,mjx,mgmyp,mkx,mgmzp,  &
       istart,iend,j,'deep',maxens2,maxens,ipr,jpr,iens,pr_ens,outt_ens,    &
       maxens3,ensdim,massfln,xfac1)

  do I=ISTART,IEND
     PRE(I) = max(PRE(I),0.)
  enddo
  !
  !---------------------done------------------------
  !
  if (j.eq.jpr) then
     i=ipr
     print*,'-----------------------done---------'
     print*,'mynum i j pre  ierr=',mynum,i,j,pre(i),ierr(i)
     do k=1,mkx
        !print *,k,z(i,k),outt(i,k),outq(i,k),outqc(i,k)
     enddo
  endif

  !---Salva parametros da CUP para uso no transporte
  !   convectivo:
  do i=istart,iend
     ierr4d(i)   = ierr(i)
     jmin4d(i)   = jmin(i)
     kdet4d(i)   = kdet(i)
     k224d(i)    = k22(i)
     kbcon4d(i)  = kbcon(i)
     ktop4d(i)   = ktop(i)
     kstabi4d(i) = kstabi(i)
     kstabm4d(i) = kstabm(i)
     !kpbl4d(i)   = kpbl(i)
     kpbl4d(i)   = k22(i) ! por enquanto kpbl==k22
     xmb4d(i)    = xmb(i)
     !print*,'1 mynum i j ierr=',mynum,i,j,ierr4d(i)
     !if(ierr(i).eq.0)  print*,'cup',ngrid,j,i,ierr4d(i),xmb4d(i)
     !srf - Media no ensemble 2 do parametro edt
     edt_average = 0.
     do iedt=1,maxens2
        edt_average = edt_average + edtc(i,iedt)
     enddo

     edt4d(i) = edt_average/float(maxens2)
     do k=1,mkx
        if (iens.eq.1) then
           zcup5d(k,i) = zo_cup(i,k)
           pcup5d(k,i) = po_cup(i,k)
        endif
        enup5d(k,i) = mentr_rate
        endn5d(k,i) = mentrd_rate
        deup5d(k,i) = cd(i,k)
        dedn5d(k,i) = cdd(i,k)
        zup5d(k,i)  = zuo(i,k)
        zdn5d(k,i)  = zdo(i,k)
        !Lufla p_lw5d(k,i) = xmb(i)*pwo(i,k)/(qrco(i,k) + 1.e-16)
        !Lufla begin -------	
	prup5d(k,i) = xmb(i)*pwo(i,k) !only for upfradt
	clwup5d(k,i) = qrco(i,k)      !only for upfradt
	tup5d(k,i) = t_cup(i,k)       !>>> em verdade deveria ser a temperatura da parcela
        !>>> de ar no updraft e _NAO_ a temperatura ambiente
        !Lufla end -------	

     enddo
  enddo

  !-------Salva parametros nas analises do RAMS

  do I=ISTART,IEND
     xierr(i,j)=float(ierr(i))
     if (ierr(i).eq.0) then
        upmf(i,j)   = xmb(i)
        !!!dnmf(i,j)   = edt4d(i)*xmb(i)
        ! downdraft mass flux averaged
	dnmf(i,j)  = 0.
        do k=1,ensdim
           dnmf(i,j) = dnmf(i,j) + massfln(i,j,k)
        enddo
        dnmf(i,j) = dnmf(i,j)/float(ensdim)
        !recalcula o parametro edt4d para o transporte convectivo 
        !modifique posteriormente para uso direto do fluxo de massa do downdraft
        !dnmf(i,j) = edt4d(i)*xmb(i)    ! downdraft mass flux averaged
        edt4d(i)  = dnmf(i,j) / ( upmf(i,j) + 1.e-16 )
 !
        xktop(i,j)  = float(ktop(i))
        xkbcon(i,j) = float(kbcon(i))
        xkdt(i,j)   = float(kdet(i))
        xjmin(i,j)  = float(jmin(i))
	xk22(i,j)   = float(k22(i))
     elseif (ierr(i).ne.0.and.ierr(i).ne.20) then
        upmf(i,j)   = 0.
        dnmf(i,j)   = 0.
        xktop(i,j)  = 0.
        xkbcon(i,j) = 0.
        xkdt(i,j)   = 0.
        xjmin(i,j)  = 0.
        xk22(i,j)   = 0.  !Lufla
     endif
  enddo

  !srf----------------------------------------

  return
end subroutine CUP_enss


!---------------------------------------------
subroutine cup_dellas(ierr, z_cup, p_cup, hcd, edt, zd, cdd, he, mix,   &
     mgmxp, mkx, mgmzp, istart, iend, della, itest, j, mentrd_rate, zu, &
     cd, hc, ktop, k22, kbcon, mentr_rate, jmin, he_cup, kdet, kpbl,    &
     ipr, jpr, name)
  implicit none
  character (LEN=*) name  !CHARACTER *(*) name
  integer mix, mgmxp, mkx, mgmzp, i, k, istart, iend,  &
       itest, j
  real z_cup(mgmxp,mgmzp), p_cup(mgmxp,mgmzp), hcd(mgmxp,mgmzp),  &
       zd(mgmxp,mgmzp), cdd(mgmxp,mgmzp), he(mgmxp,mgmzp),        &
       della(mgmxp,mgmzp), hc(mgmxp,mgmzp), cd(mgmxp,mgmzp),      &
       zu(mgmxp,mgmzp), he_cup(mgmxp,mgmzp)
  real edt(mgmxp)
  integer kbcon(mgmxp), ktop(mgmxp), k22(mgmxp), jmin(mgmxp)
  integer ierr(mgmxp), kdet(mgmxp), kpbl(mgmxp), ipr, jpr
  real detdo1, detdo2, entdo, g, dp, dz, mentrd_rate,  &
       mentr_rate, subin, detdo, entup, detup,         &
       subdown, entdoj, entupk, detupk, totmas
  real xsum, xsumt
  g = 9.81
  i = istart
  if (j.eq.jpr) then
     if (itest.eq.1) then
        print *, 'this one for h '
        print *, 'in dellas kpbl(i),k22(i),kbcon(i),ktop(i),jmin(i)'
        print *, kpbl(i), k22(i), kbcon(i), ktop(i), jmin(i)
     else
        print *, 'this one for q '
        print *, 'in dellas kpbl(i),k22(i),kbcon(i),ktop(i),jmin(i)'
        print *, kpbl(i), k22(i), kbcon(i), ktop(i), jmin(i)
     endif
  endif
  do K=2,MKX
     do I=ISTART,IEND
        della(i,k) = 0.
     enddo
  enddo
  xsum=0.
  !      hesum=0.
  xsumt=0.

  do K=2,MKX-1
  !DO 100 K=2,MKX-1
     do I=ISTART,IEND
     !DO 100 I=ISTART,IEND
        !IF (ierr(i).NE.0) GO TO 100
        if (ierr(i).ne.0) cycle
        !IF (K.GT.KTOP(I)) GO TO 100
        if (K.gt.KTOP(I)) cycle
        !
        !--- Specify detrainment of downdraft,
        !    has to be consistent
        !--- with zd calculations in soundd.
        !
        dz    = Z_cup(I,K+1)-Z_cup(I,K)
        detdo = edt(i)*CDD(i,K)   *dz*zd(i,k+1)
        entdo = edt(i)*mentrd_rate*dz*zd(i,k+1)
        subin = zu(i,k+1)-zd(i,k+1)*edt(i)
        entup = 0.
        detup = 0.
        if (k.ge.kbcon(i).and.k.lt.ktop(i)) then
           entup = mentr_rate*dz*zu(i,k)
           detup = CD(i,K+1) *dz*zu(i,k)
        endif
        subdown = ( zu(i,k)-zd(i,k)*edt(i) )
        entdoj  = 0.
        entupk  = 0.
        detupk  = 0.

        if (k.eq.jmin(i)) then
           entdoj  = zd(i,k)*edt(i)
        endif

        if (k.eq.k22(i)-1) then
           !        IF (k.EQ.kpbl(i)) THEN
           entupk  = zu(i,kpbl(i))
        endif

        if (k.gt.kdet(i)) then
           detdo   = 0.
        endif

        if (k.eq.ktop(i)-0) then
           detupk  = zu(i,ktop(i))
           subin   = 0.
        endif
        if (k.lt.kbcon(i)) then
           detup   = 0.
        endif
        !
        !--- Changed due to subsidence and entrainment
        !
        totmas =subin-subdown+detup-entup-entdo + detdo-entupk-entdoj+detupk
        if (j.eq.jpr.and.i.eq.ipr)                       &
             print *, 'k,totmas,sui,sud = ', k, totmas, subin,subdown
        if (j.eq.jpr.and.i.eq.ipr)                       &
             print *, 'updr stuff = ', detup, entup, entupk, detupk
        if (j.eq.jpr.and.i.eq.ipr)                       &
             print *, 'dddr stuff = ', entdo, detdo, entdoj
        if (abs(totmas).gt.1.e-6) then
           print *, '*********************', i, j, k, totmas, name
           print *, kpbl(i), k22(i), kbcon(i), ktop(i)
           !          print *,'updr stuff = ',subin,
           !    1      subdown,detup,entup,entupk,detupk
           !          print *,'dddr stuff = ',entdo,
           !    1      detdo,entdoj
           stop
        endif

        !srf         dp =  100.*( p_cup(i,k-1)-p_cup(i,k) )
        dp =  100.*( p_cup(i,k)-p_cup(i,k+1) )
        della(i,k)=(subin  *he_cup(i,k+1) - subdown*he_cup(i,k  ) +  &
             detup*.5*( HC(i,K+1)+ HC(i,K)) +                        &
             detdo*.5*(HCD(i,K+1)+HCD(i,K)) -                        &
             entup*he(i,k) - entdo*he(i,k) -                         &
             entupk*he_cup(i,k22(i)) -  entdoj*he_cup(i,jmin(i)) +   &
             detupk*hc(i,ktop(i)))*g/dp

        !     if(j.eq.3.and.i.eq.120)xsumt=xsumt+totmas
        !     if(j.eq.3.and.i.eq.120)hesum=hesum+he(i,k)*dp
        !     if(j.eq.3.and.i.eq.120)xsum=xsum+della(i,k)*dp
        !     if(j.eq.3.and.i.eq.120)print *,'xsum = ',xsum
        if (i.eq.ipr.and.j.eq.jpr) then
           print *, k, della(i,k), subin*he_cup(i,k+1),                  &
                subdown*he_cup(i,k), detdo*.5*(HCD(i,K+1)+HCD(i,K))
              print *, k, detup*.5*(HC(i,K+1)+HC(i,K)),                  &
                   detupk*hc(i,ktop(i)), entup*he(i,k),entdo*he(i,k)
              print *, k, he_cup(i,k+1), he_cup(i,k), entupk*he_cup(i,k)
           endif

!100     CONTINUE
     enddo
  enddo

  return
end subroutine cup_dellas


!-----------------------------------------
subroutine cup_dellabot(ipr, jpr, he_cup, ierr, z_cup, p_cup, hcd, edt,   &
     zd, cdd, he, mix, mgmxp, mkx, mgmzp, istart, iend, della, itest, j,  &
     mentrd_rate, z)
  implicit none
  integer mix, mgmxp, mkx, mgmzp, i, istart, iend, itest, j
  real z_cup(mgmxp,mgmzp), p_cup(mgmxp,mgmzp), hcd(mgmxp,mgmzp),   &
       zd(mgmxp,mgmzp), cdd(mgmxp,mgmzp), he(mgmxp,mgmzp),         &
       della(mgmxp,mgmzp), he_cup(mgmxp,mgmzp), z(mgmxp,mgmzp), edt(mgmxp)
  integer ierr(mgmxp), ipr, jpr, m
  real detdo1, detdo2, entdo, g, dp, dz, mentrd_rate, subin, detdo

  g = 9.81
  do i=istart,iend
  !DO 100 i=istart,iend
     della(i,1)=0.
     !IF (ierr(i).NE.0) GO TO 100
     if (ierr(i).ne.0) cycle
     dz        =       z_cup(i,2)-z_cup(i,1)
     dp        = 100.*(p_cup(i,1)-p_cup(i,2))
     detdo1    = edt(i)*zd(i,2)*cdd(i,1)*dz
     detdo2    = edt(i)*zd(i,1)
     entdo     = edt(i)*zd(i,2)*mentrd_rate*dz
     subin     =-edt(I)*zd(i,2)
     detdo     = detdo1+detdo2-entdo+subin
     DELLA(I,1)= (detdo1*.5*(hcd(i,1)+hcd(i,2)) +  &
          detdo2*hcd(i,1) + subin*he_cup(i,2) -    &
          entdo*he(i,1))*g/dp
     !
     !if(i.eq.ipr.and.j.eq.jpr)print*,detdo1,detdo2,entdo,
     !1     subin,detdo
     !if(i.eq.ipr.and.j.eq.jpr)print *,della(i,1),hcd(i,1),hcd(i,2),
     !1     he_cup(i,2),he(i,1)
     !srf-----print-------
     if (j.eq.jpr   .and. i.eq.ipr) then
        write(6,'(a1,78a1)') ' ',('-',m=1,78)
        print*, 'i k Z_cup1 P_cup1 Z_cup2 P_cup2 dp dz dens della '
        write(6,'(2i4,8F12.4)') i, 1, Z_cup(i,1),                             &
             p_cup(i,1), Z_cup(i,2), p_cup(i,2),                              &
             100.*(p_cup(i,1)-p_cup(i,2)), dz, 100.*(p_cup(i,1)-p_cup(i,2))/  &
             (g*(z_cup(i,2)-z_cup(i,1))), della(i,1)
        ! densidade

     endif
     !srf-----print-------

!100  CONTINUE
  enddo
  return
end subroutine cup_dellabot


!------------------------------------------------------------
subroutine cup_forcing_ens_16(aa0,aa1,xaa0,mbdt,dtime,xmb,ierr,   &
     mix,mgmxp,mjx,mgmyp,mkx,mgmzp,istart,iend,xf,j,fquasi,       &
     fstab,name,xland,maxens,iens,iedt,maxens2,ipr,jpr,maxens3,   &
     mconv,omeg,zd,k22,zu,pr_ens,edt,aad,kbcon,massflx,		  &
     iact_old_gr,dir,ensdim,massfln,massfld,iresult,xff_ens3,xk,  &
     p_cup,ktop,icoic)

  implicit none
  character (LEN=*) name  !CHARACTER *(*) name

  integer k,i,istart,iend,mix,mgmxp,mjx,mgmyp,mkx,mgmzp,j,        &
       maxens,maxens3
  integer ensdim,iens,nall,iedt,maxens2,ipr,jpr

  !------ ensemble 3 dimension = 16
  integer mkxcrt,kclim
  parameter (mkxcrt=15)
  real pcrit(mkxcrt),acrit(mkxcrt),acritt(mkxcrt),aclim1,         &
       aclim2,aclim3,aclim4
  data pcrit/850.,800.,750.,700.,650.,600.,550.,500.,450.,400.,   &
       350.,300.,250.,200.,150./
  data acrit/.0633,.0445,.0553,.0664,.075,.1082,.1521,.2216,      &
       .3151,.3677,.41,.5255,.7663,1.1686,1.6851/
  !  GDAS derived acrit
  data acritt/.203,.515,.521,.566,.625,.665,.659,.688,            &
       .743,.813,.886,.947,1.138,1.377,1.896/
  integer ktop(mgmxp)
  real p_cup(mgmxp,mgmzp)
  !------

  integer k22(mgmxp),kbcon(mgmxp),ierr(mgmxp)
  integer iact_old_gr(mgmxp,mgmyp)

  real aa0(mgmxp),aa1(mgmxp),xaa0(mgmxp,maxens),xmb(mgmxp)
  real mbdt(maxens),dtime,dxxf,edt(mgmxp),aad(mgmxp),dir(mgmxp)
  real     xf(mgmxp,mgmyp,ensdim),xland(mgmxp,mgmyp)
  real pr_ens(mgmxp,mgmyp,ensdim)
  real mconv(mgmxp),omeg(mgmxp,mgmzp),zd(mgmxp,mgmzp),            &
       zu(mgmxp,mgmzp)
  real xff_ens3(maxens3),xk(maxens),xff,xff1,xff2,xff3,xff0
  real massflx(mgmxp,mgmyp)
  real massfln(mgmxp,mgmyp,ensdim)
  real xomg,massfld
  integer fquasi,fstab,nens,ne,n,nens3,iresult,iresultd,          &
       iresulte,icoic
  nens=0

  !--- LARGE SCALE FORCING
  !
  !DO 100 I=ISTART,IEND
  do I=ISTART,IEND
     xmb(i)=0.
     if(name.eq.'deeps'.and.ierr(i).gt.995)then
        aa0(i) =0.
        ierr(i)=0
     endif
     if(ierr(i).eq.0)then
        !Added for ensemble 3 with dimension = 16
        kclim=0
        do k=mkxcrt,1,-1
           if(p_cup(i,ktop(i)).lt.pcrit(k))then
              kclim=k
              go to 9
           endif
        enddo
        if(p_cup(i,ktop(i)).gt.pcrit(1))kclim=1
9       continue
        k= max(kclim-1,1)
        aclim1= acrit(kclim)*1.e3
        aclim2= acrit(k)*1.e3
        aclim3= acritt(kclim)*1.e3
        aclim4= acritt(k)*1.e3
        !
        !--- Treatment different for this closure
        !
        if(name.eq.'deeps')then
           !
           xff0       =  (AA1(I)-AA0(I))/dtime
           xff_ens3(1)=  (AA1(I)-AA0(I))/dtime
           xff_ens3(2)= .9*xff_ens3(1)
           xff_ens3(3)=1.1*xff_ens3(1)
           !
           !     More like Brown (1979), or Frank-Cohen (199?)
           !
           !---  omeg is in bar/s, mconv done with omeg in Pa/s
           xff_ens3(4)=     -omeg(i,k22(i))/9.81
           xff_ens3(5)=     -omeg(i,kbcon(i))/9.81
           xff_ens3(6)=     -omeg(i,1)/9.81
           do k=2,kbcon(i)-1
              xomg=     -omeg(i,k)/9.81
              if(xomg.gt.xff_ens3(6)) xff_ens3(6)=xomg
           enddo
           !
           !--- More like Krishnamurti et al.
           !
           xff_ens3(7)=    mconv(i)
           xff_ens3(8)= .9*mconv(i)
           xff_ens3(9)=1.1*mconv(i)
           !
           !--- More like Fritsch Chappel or Kain Fritsch (plus triggers)
           !srf - changed at dec/2002 - greater timescale instab. removal
           !xff_ens3(10)=AA1(I)/(60.*20.)
           !xff_ens3(11)=AA1(I)/(60.*30.)
           !xff_ens3(12)=AA1(I)/(60.*40.)
           xff_ens3(10)=AA1(I)/(60.*50.)
           xff_ens3(11)=AA1(I)/(60.*60.)
           xff_ens3(12)=AA1(I)/(60.*70.)

           !   
           !--- More original Arakawa-Schubert (climatologic value of aa0)
           !
           xff_ens3(13)=max(0.,(AA1(I)-aclim1)/dtime)
           xff_ens3(14)=max(0.,(AA1(I)-aclim2)/dtime)
           xff_ens3(15)=max(0.,(AA1(I)-aclim3)/dtime)
           xff_ens3(16)=max(0.,(AA1(I)-aclim4)/dtime)

           do nens=1,maxens
              XK(nens)=(XAA0(I,nens)-AA1(I))/MBDT(nens)
              if(xk(nens).le.0.and.xk(nens).gt.-1.e-9) xk(nens)=-1.e-9
              if(xk(nens).gt.0.and.xk(nens).lt.+1.e-9) xk(nens)=+1.e-9
           enddo
           !
           !--- Add up all ensembles
           !
           do ne=1,maxens
              !
              !--- for every xk, we have maxens3 xffs
              !--- iens is from outermost ensemble (most expensive!
              !
              !--- iedt (maxens2 belongs to it)
              !--- is from second, next outermost, not so expensive
              !
              !--- so, for every outermost loop, we have maxens*maxens2*3
              !--- ensembles!!! nall would be 0, if everything is on first
              !--- loop index, then ne would start counting, then iedt, 
              !--- then iens....
              !
              iresultd=0
              iresulte=0
              nall=(iens-1)*maxens3*maxens*maxens2  &
                   +(iedt-1)*maxens*maxens3         &
                   +(ne-1)*maxens3
              !
              !--- check for upwind convection
              !
              !iresult=0
              !massfld=0.
              !call cup_direction2(i,j,dir,iact_old_gr,mix,mjx,  &
              !     mgmxp,mgmyp,massflx,iresult,ensdim,1,nall,   &
              !     maxens3,massfld)
              !if(i.eq.ipr.and.j.eq.jpr.and.iedt.eq.1.and.ne.eq.1)then
              !   print *,massfld,ne,iedt,iens
              !   print *,xk(ne),xff_ens3(1),xff_ens3(2),xff_ens3(3)
              !endif
              if(XK(ne).lt.0.and.xff0.gt.0.)iresultd=1
              iresulte=max(iresult,iresultd)
              iresulte=1
              if(iresulte.eq.1)then
                 !
                 !--- Special treatment for stability closures
                 !

                 if(xff0.gt.0.)then
                    xf(i,j,nall+1) =max(0., -xff_ens3(1)/xk(ne))+massfld
                    xf(i,j,nall+2) =max(0., -xff_ens3(2)/xk(ne))+massfld
                    xf(i,j,nall+3) =max(0., -xff_ens3(3)/xk(ne))+massfld
                    xf(i,j,nall+13)=max(0.,-xff_ens3(13)/xk(ne))+massfld
                    xf(i,j,nall+14)=max(0.,-xff_ens3(14)/xk(ne))+massfld
                    xf(i,j,nall+15)=max(0.,-xff_ens3(15)/xk(ne))+massfld
                    xf(i,j,nall+16)=max(0.,-xff_ens3(16)/xk(ne))+massfld
                 else
                    xf(i,j,nall+1) =massfld
                    xf(i,j,nall+2) =massfld
                    xf(i,j,nall+3) =massfld
                    xf(i,j,nall+13)=massfld
                    xf(i,j,nall+14)=massfld
                    xf(i,j,nall+15)=massfld
                    xf(i,j,nall+16)=massfld
                 endif
                 !
                 !--- if iresult.eq.1, following independent of xff0
                 !
                 xf(i,j,nall+4)=max(0.,xff_ens3(4)+massfld)
                 xf(i,j,nall+5)=max(0.,xff_ens3(5)+massfld)
                 xf(i,j,nall+6)=max(0.,xff_ens3(6)+massfld)
                 xf(i,j,nall+7)=max(0.,xff_ens3(7)/pr_ens(i,j,nall+7))
                 xf(i,j,nall+8)=max(0.,xff_ens3(8)/pr_ens(i,j,nall+8))
                 xf(i,j,nall+9)=max(0.,xff_ens3(9)/pr_ens(i,j,nall+9))
                 if(XK(ne).lt.0.)then
                    xf(i,j,nall+10)=max(0.,-xff_ens3(10)/xk(ne))+massfld
                    xf(i,j,nall+11)=max(0.,-xff_ens3(11)/xk(ne))+massfld
                    xf(i,j,nall+12)=max(0.,-xff_ens3(12)/xk(ne))+massfld
                 else
                    xf(i,j,nall+10)=massfld
                    xf(i,j,nall+11)=massfld
                    xf(i,j,nall+12)=massfld
                 endif
                 if(icoic.ge.1)then
                    xf(i,j,nall+1) =xf(i,j,nall+icoic)
                    xf(i,j,nall+2) =xf(i,j,nall+icoic)
                    xf(i,j,nall+3) =xf(i,j,nall+icoic)
                    xf(i,j,nall+4) =xf(i,j,nall+icoic)
                    xf(i,j,nall+5) =xf(i,j,nall+icoic)
                    xf(i,j,nall+6) =xf(i,j,nall+icoic)
                    xf(i,j,nall+7) =xf(i,j,nall+icoic)
                    xf(i,j,nall+8) =xf(i,j,nall+icoic)
                    xf(i,j,nall+9) =xf(i,j,nall+icoic)
                    xf(i,j,nall+10)=xf(i,j,nall+icoic)
                    xf(i,j,nall+11)=xf(i,j,nall+icoic)
                    xf(i,j,nall+12)=xf(i,j,nall+icoic)
                    xf(i,j,nall+13)=xf(i,j,nall+icoic)
                    xf(i,j,nall+14)=xf(i,j,nall+icoic)
                    xf(i,j,nall+15)=xf(i,j,nall+icoic)
                    xf(i,j,nall+16)=xf(i,j,nall+icoic)
                 endif
                 !==============
                 !05-12-2002
                 !srf - forcing 14 is too bad, use the same for 13:
                 !A&S (14) = A&S (13)
                 xf(i,j,nall+14)=xf(i,j,nall+13)
                 !==============
                 !
                 !--- store new for next time step
                 !
                 do nens3=1,maxens3
                    massfln(i,j,nall+nens3)=edt(i)*xf(i,j,nall+nens3)
                    massfln(i,j,nall+nens3)=max(0.,massfln(i,j,nall+nens3))
                 enddo
              endif
           enddo
           cycle !go to 100
        endif
     elseif(ierr(i).ne.20.and.ierr(i).ne.0)then
        do n=1,ensdim
           xf(i,j,n)=0.
           massfln(i,j,n)=0.
        enddo
     endif
     !100  CONTINUE
  enddo
  return
end subroutine cup_forcing_ens_16
!--------------------------------------------------------------------

!--------------------------------------------------------------------
subroutine cup_output_ens(xf,ierr,dellat,dellaq,dellaqc,    &
     outtem,outq,outqc,pre,pw,xmb,ktop,mix,mgmxp,mjx,mgmyp,   &
     mkx,mgmzp,istart,iend,j,name,nx,nx2,ipr,jpr,iens,pr_ens, &
     outt_ens,maxens3,ensdim,massfln,xfac1)
  implicit none
  character (LEN=*) name  !CHARACTER *(*) name

  integer mix,mjx,mkx,istart,iend,mgmxp,mgmyp,mgmzp, &
       ensdim,i,k,j,nx,n,nx2,ipr,jpr,m

  real       xf(mgmxp,mgmyp,ensdim), pr_ens(mgmxp,mgmyp,ensdim), &
       outt_ens(mgmxp,mgmyp,ensdim),massfln(mgmxp,mgmyp,ensdim)

  real outtem(mgmxp,mgmzp),outq(mgmxp,mgmzp),outqc(mgmxp,mgmzp), &
       dellat(mgmxp,mgmzp,nx), dellaq(mgmxp,mgmzp,nx),           &
       pw(mgmxp,mgmzp,nx),dellaqc(mgmxp,mgmzp,nx),               &
       pre(mgmxp),xmb(mgmxp),xfac1(mgmxp)

  integer ktop(mgmxp),ierr(mgmxp),ncount,iens,maxens3,nens3
  real outtes,ddtes,dtt,dtq,dtqc,dtpw

  do K=1,MKX
     do I=ISTART,IEND
        outtem(i,k) = 0.
        outq(i,k) = 0.
        outqc(i,k) = 0.
     enddo
  enddo

  do I=ISTART,IEND
     pre(i)  =0.
     xmb(i)  =0.
     xfac1(i)=1.
  enddo
  !
  !--- Calculate ensemble average mass fluxes
  !
  do I=ISTART,IEND
     ncount=0
     xmb(i)=0.

     if(ierr(i).eq.0)then

        do n=(iens-1)*nx*nx2*maxens3+1,iens*nx*nx2*maxens3
           pr_ens(i,j,n) =   pr_ens(i,j,n)*xf(i,j,n)
           outt_ens(i,j,n) = outt_ens(i,j,n)*xf(i,j,n)

           if(xf(i,j,n).gt.0.)then
              xmb(i) = xmb(i) + xf(i,j,n)
              ncount = ncount + 1 
              !if(i.eq.ipr.and.j.eq.jpr) print *,'XF =',n,xf(i,j,n)
           endif
        enddo

        if(ncount.gt.0)then
           xmb(i)=xmb(i)/float(ncount)
        else
           xmb(i)=0.
           ierr(i)=13
        endif

     endif

  enddo
  !
  !-- Now do feedback
  !
  ddtes=250.
  !     if(name.eq.'shal')ddtes=500.
  do I=ISTART,IEND
     xfac1(i)=xmb(i)
  enddo
  do K=1,MKX
     do I=ISTART,IEND
        dtt  =0.
        dtq  =0.
        dtqc =0.
        dtpw =0.

        if(ierr(i).eq.0.and.k.le.ktop(i))then
           do n=1,nx
              dtt  = dtt  +  dellat(i,k,n)
              dtq  = dtq  +  dellaq(i,k,n)
              dtqc = dtqc + dellaqc(i,k,n)
              dtpw = dtpw +      pw(i,k,n)

              !srf-----print-------
              !	   if(k.eq.1) then
              !            write(6,'(a1,78a1)') ' ',('*',m=1,78)
              !	    print*,'nx j i k PREC dtpw  pw '
              !	   endif
              !          write(6,'(4i4,3e12.4)') nx,j,i,k,pre(i),dtpw,pw(i,k,n)
              !     	   if(k.eq.ktop(i)) then
              !            write(6,'(a1,78a1)') ' ',('*',m=1,78)
              !           endif
              !srf-----print-------
           enddo
           outtes = dtt*XMB(I)*86400./float(nx)

           if (outtes .gt. 2.*ddtes .and. k.gt.2) then
              XMB(I) = 2.*ddtes/outtes * xmb(i)
              outtes = 1.*ddtes
           endif

           if (outtes .lt. -ddtes)                then
              XMB(I) = -ddtes/outtes * xmb(i)
              outtes = -ddtes
           endif

           if (outtes .gt. .5*ddtes .and. k.le.2) then
              XMB(I) =    ddtes/outtes * xmb(i)
              outtes = .5*ddtes
           endif

           OUTTEM(I,K) = OUTTEM(I,K) +XMB(I)*dtt /float(nx)
           OUTQ(I,K) =   OUTQ(I,K) +XMB(I)*dtq /float(nx)
           OUTQC(I,K) =  OUTQC(I,K) +XMB(I)*dtqc/float(nx)
           PRE(I)      = PRE(I)      +XMB(I)*dtpw/float(nx) 
           ! unit : kg[liq water]/(m^2 s)

           !srf-----print-------
           !          if(j.eq.jpr   .and. i.eq.ipr) then
           !	   if(k.eq.1) then
           !            write(6,'(a1,78a1)') ' ',('-',m=1,78)
           !	    print*,'nx j i k OUTTEMP  OUTQ OUTQC PREC dtpw xmb pw'
           !	   endif
           !          write(6,'(4i4,8e12.4)') nx,j,i,k
           !     &   ,86400.*OUTTEM(I,K),1000.*86400.*OUTQ(I,K)
           !     &   ,1000.*86400.*OUTQC(I,K)
           !     &   ,pre(i),dtpw,xmb(i),pw(i,k,n)
           !     
           !
           !	   if(k.eq.ktop(i)) then
           !            write(6,'(a1,78a1)') ' ',('-',m=1,78)
           !           endif
           !           if(k.eq.mkx) write(6,'(a1,78a1)') ' ',('-',m=1,78)
           !          endif
           !srf-----print-------

        endif
     enddo
  enddo

  do I=ISTART,IEND
     if(ierr(i).eq.0)then
        xfac1(i)=xmb(i)/xfac1(i)
        do k=1,ensdim
           massfln(i,j,k)=massfln(i,j,k)*xfac1(i)
        enddo
     endif
  enddo

  return
end subroutine cup_output_ens
!------------------------------------------------------------
